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The three T states, T(1S+2S+3S), are measured in d+Au and p+p collisions at ^/s NN —200 GeV 
and rapidities 1.2 < \y\ < 2.2 by the PHENIX experiment at the Relativistic Heavy-Ion Collider. 
Cross sections for the inclusive Y(1S+2S+3<S') production are obtained. The inclusive yields per 
binary collision for d+Au collisions relative to those in p+p collisions (RdAu) are found to be 0.62 
± 0.26 (stat) ± 0.13 (syst) in the gold-going direction and 0.91 ± 0.33 (stat) ± 0.16 (syst) in the 
deuteron-going direction. The measured results are compared to a nuclear-shadowing model, EPS09 
[JHEP 04, 065 (2009)], combined with a final-state breakup cross section, at r , and compared to 
lower energy p+A results. We also compare the results to the PHENIX J/ip results [Phys. Rev. 
Lett. 107, 142301 (2011)]. The rapidity dependence of the observed T suppression is consistent 
with lower energy p+A measurements. 

PACS numbers: 25.75.Dw 



I. INTRODUCTION 

Quarkonia are produced dominantly by the gluon- 
gluon fusion process in high energy collisions Jl|, 
Therefore, quarkonia production is a good probe to ex- 
plore the gluon distribution of the nucleon and its modifi- 
cation in nuclei. Recently, the PHENIX collaboration has 
reported J/ip suppression in ^/s NN —200 GeV deuteron- 
gold (d+Au) collisions at the Relativistic Heavy Ion Col- 
lider 3]. The centrality dependence of these J/ip sup- 
pression results at forward rapidity is not well described 
quantitatively by nuclear-shadowing models that include 
final-state breakup effects [4]. Because the T mass is 
heavier than J/ip, the nuclear effects on the gluon distri- 
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bution can be studied in different kinematic regions. At 
forward rapidity (the deuteron going direction) and the 
same collision energy of y/s^ =200 GeV, the average 
momentum fraction of the gluon in the gold nucleus that 
is sampled for T production is (x 2 ) ~lx 10~ 2 , whereas 
J/ip production samples (x 2 ) ~ 3 x 10~ 3 . 

There are various fits for the nuclear parton distribu- 
tion functions (nPDFs) over broad x ranges In 
d+Au collisions, since the forward and backward rapidi- 
ties cover different x 2 (x in the Au nucleus) ranges, T 
production at these two rapidities would be affected dif- 
ferently by these nPDFs. Additionally, the final-state 
breakup effect should also suppress T yields by some 
amount at both rapidities, but there is no clear indication 
of the size of this effect yet Thus, T measurements 
in <i+Au collisions should give new and valuable infor- 
mation to test nuclear parton modification and breakup 
effects. 

Lattice quantum chromodynamics predicts that the 
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T(15), T(2S), and T(35) all have different binding en- 
ergies and radii, and so should melt at different tem- 
peratures of the hot nuclear medium (Tlj . Therefore, 
the three T states are thought to be good probes for 
the temperature of the hot dense matter. Recently, the 
Compact Muon Solenoid (CMS) experiment at the Large 
Hadron Collider (LHC) reported that the double ratio of 
the T(2S+3S) excited states to the T(1S) ground state 
in Pb+Pb and p+p collisions at y/s NN —2.76 TeV, 



2008 



T(2S+3S)/T(lS)| pb+pb 
T(2S+3S)/T(lS)\ p+p 



= 0.31 



+0.19 
-0.15 



(stat) ± 0.03(syst) 



for single decay muons of px > 4 GeV/c and \t)\ < 2.4 [l^ . 
They also reported T(1S) is suppressed by approximately 
40% in minimum-bias Pb+Pb collisions (12| . 

In T suppression for nucleus-nucleus collisions, there 
should be contributions from cold nuclear matter as well 
as those from the hot nuclear matter. Thus, to separate 
these two types of contributions, it is necessary to mea- 
sure the level of suppression from cold nuclear matter 
effects with p(d)+A collisions, where hot nuclear matter 
is not created. 

A lower-energy fixed-target experiment, E772, re- 
ported measurements in y/s NN =38.8 GeV p + A col- 
lisions of the T(15) ground state and the T(2S+35) 
excited states. The observed suppression of the T(15) 
and T(2S+3S') agree within the experimental uncertain- 
ties [13[ • The initial-state effects from nuclear shadowing 
are not expected to differ between the three T states since 
they are produced mostly by gluon-gluon fusion sub-pro- 
cesses and have similar masses [l|, 0, [l(| • For the final- 
state breakup effect, there is no clear estimate of its en- 
ergy dependence and of the difference between the three 
T states. In this paper, we present the first measurement 
of inclusive T(1S+2S+3S) cold nuclear matter effects as 
well as the production cross section using <i±Au and p+p 
collisions at y/J^ —200 GeV measured by the PHENIX 
experiment. 



II. ANALYSIS METHOD 
A. Experimental setup 

The PHENIX apparatus is described in detail in 
Ref. In d+Au collisions, the deuteron comes from 

the negative-rapidity end of PHENIX (South) and goes 
towards positive-rapidity (North) , and vice versa for the 
gold ions. For the T analysis presented here, three de- 
tector systems are required for reconstruction and trig- 
gering at forward and backward rapidities. These are the 
Muon Tracker (MuTr), the Muon Identifier (MuID), and 
the beam-beam counters (BBCs). There are two sepa- 
rate BBC systems. One covers forward rapidity and the 
other covers backward rapidity as shown in Fig. [1] 

Each BBC comprises 64 quartz Cerenkov radiators 
and mesh dynode PMTs. The two BBCs are located at 
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FIG. 1: (color online) PHENIX detector configuration in 
2008. This side view includes the forward-rapidity detec- 
tors (South and North arms): the Muon Tracker (MuTr), 
the Muon Identifier (MuID) for muon detection and identifi- 
cation, and the beam-beam counter (BBC) for global event 
characteristics. 



±144 cm from the nominal interaction point and cover 
pseudorapidity of 3 < \rj\ < 3.9. Each MuID comprises 
five layers of thick steel vertical plates with Iarocci tubes 
between each pair of plates. Most hadrons are absorbed 
in the steel plates. Muons with more than 2.7 GeV/c 
of momentum will pass through all layers of the MuID 
and reach the last gap. Each MuTr is composed of three 
stations of cathode strip chambers and measures the mo- 
mentum and charge sign of the muon according to their 
bending in the magnetic field, with coverage in rapidity 
of 1.2 < \y\ < 2.2 for the T and full azimuthal coverage 
of <f> £ [— 7T, tt]. The nose-cone absorber and the central- 
magnet pole face, which both lie between the interaction 
region and the innermost part of the muon tracker, also 
help to reduce hadron backgrounds, especially by elimi- 
nating many light hadrons (e.g. tt, K) before they decay 
into secondary muons. Fewer than 1% of hadrons punch 
through the absorbers, reach the last gap of MuID, and 
become fake muon tracks. 

The data sets used in this analysis were collected dur- 
ing 2006, 2008, and 2009 using the BBC Level- 1 trig- 
ger. This BBC trigger requires hits in the negative and 
positive rapidity ends of the BBC in order to register 
an interaction and provide a minimum-bias trigger. The 
BBC also measures the z-vertex position of the interac- 
tion using time differences between its hits in the negative 
and positive rapidity directions. For this analysis, the z- 
vertex is required to be within ±30 cm of the center of 
PHENIX, z = 0. Additionally, the MuID Level-1 trig- 
ger is used in order to require that at least two particles 
penetrate through the MuID to its last layer. 

After removing bad runs, such as those with numerous 
high-voltage trips and significant detector performance 
variations, the integrated luminosities of the collected 
data are 69 nb _1 and 67 nb _1 for the positive and neg- 
ative rapidity muon detectors in d+Au collisions from 
2008. Here, 69 nb _1 and 67 nb _1 correspond to 27.2 pb _1 
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and 26.4 pb _1 when scaled by the number of partici- 
pants. For the p+p collisions, the integrated luminosities 
are 22.5 pb _1 and 22.2 pb _1 for the positive and negative 
rapidities from 2006 and 2009. 

We apply quality-assurance cuts on the data to select 
good tracks and improve the signal-to-background ratio. 
We calculate the track y 2 an d vertex Xvtx> an d match 
the tracks in the MuID and the MuTr at the first layer of 
the MuID in both position and angle. We also check the 
number of hits in a MuID road, which is a straight line 
that connects sets of hits in different layers of the MuID. 
We compare momenta of the two muons and remove pairs 
with a large asymmetry (\{pi — Pzj/ipi + P2)| > 0.6) 
between the two momenta. These asymmetric-momenta 
pairs are largely from random pairs where one hadron has 
decayed into a muon inside the tracking volume and has 
been misreconstructed as a higher momentum track; thus 
yielding a fake high-mass pair. The efficiency loss from 
this cut for Ts is less than 2%. The values of the cuts 
are determined using the PHENIX GEANT3-based [H| 
(PISA) detector simulations. 

For this analysis, we form an invariant-mass distri- 
bution from the unlike-charge-sign (foreground) pairs of 
muon tracks. In addition to the quarkonia resonances in- 
cluding the T signal, the mass distribution also contains 
uncorrelated (combinatorial) background and correlated 
background pairs. There are two methods to estimate 
the combinatorial backgrounds: 1) use like-sign pairs of 
muons from the same event, or 2) use an event-mixing 
method which mixes unlike-sign muons from different 
events to form random pairs. In this analysis, we use 
the event-mixing method to estimate the combinatorial 
background as shown in Fig. [5] (a), (c) and (e), and assign 
a systematic uncertainty based on the difference between 
the two methods. We calculate the normalization factor 
for the mixed events by 



Normalization Factor 



2 x yjFG++ x FG— 
BG+- ' 



(1) 



Here, BG^ stands for the number of the unlike-sign 

mixed events, FG++ and FG represent the number of 

the like-sign events. Unlike-sign mixed events are scaled 
by the normalization factor and we assign a 3% system- 
atic uncertainty for this factor. 

After the combinatorial background is subtracted, in 
the T-mass region, there are still contributions from 
correlated backgrounds expected from the Drell-Yan 
process and pairs of muons from the same cc or bb pairs. 
Therefore, it is important to estimate the correlated 
backgrounds properly to extract the T signal. We use 
next-to-leading order (NLO) calculations and pythia 
6.4 [l6[ to estimate these correlated backgrounds, and 
the PHENIX GEANT3 simulation to include a realistic 
detector response as shown in Fig. [2] (b), (d) and (f). 
Details of the estimates for these correlated backgrounds 
are described in the following sections. 



B. T and physical background estimation 

1. The Drell-Yan process 

The mass region between 4 and 8 GeV/c 2 (above the 
J/tp, tp masses and below the T mass) is dominated 
by the Drell-Yan process and by correlated open-heavy 
flavor pairs. The very low statistics above the T mass, 
where the Drell-Yan process dominates, does not provide 
a useful constraint on the Drell-Yan yield; so, we use 
NLO calculations from Vitev 17[ to constrain the Drell- 
Yan yields and to estimate their contribution in the T- 
mass region. NLO calculations of the Drell-Yan process 
are known to be very accurate from comparisons to data 
at other energies [13, [H, H3| ■ For d+Au collisions, nu- 
clear effects are added in the NLO calculations - including 
isospin effects which account for the composition of the 
nucleus in terms of neutrons and protons, parton shad- 
owing corrections, and the effect of initial-state energy 
loss [17j . 

To evaluate the model's systematic uncertainty for 
the Drell-Yan contribution, we use a calculation from 
CTEQ Jp, as shown in Fig. [3k.. For the Vitev calcu- 
lation without nuclear corrections, the difference be- 
tween the CTEQ calculation and that from Vitev is ap- 
proximately 10% over the entire mass range (filled-green 
squares and the black circles in Fig. [3_p). We assume 
this same systematic uncertainty for p+p and d+Au col- 
lisions. Additionally, the variation of the renormalization 
and factorization scales in the calculation from Vitev, 
Q/2 < /.i < 2Q, is included as a systematic uncertainty 
for both collision systems. 

The Drell-Yan contribution to the data is determined 
using the calculated cross section and the integrated lu- 
minosity for each data set. This contribution is corrected 
for geometrical acceptance and efficiencies. Details are 
shown in Eq. El and dU; and are discussed in the text 
that follows. 



deny r bbc 



Ni 



MB 



dar>Y 



dm cr Tot • e^BC 



dauY N M b 



e DY C ' ^ e DY — 



V, 



DY 



Am ' 



dm OTot • C 



AevY 



Nby 
Am ' 



(2) 



(3) 



(4) 



where c?ctdy /dm is the differential cross section of the 
Drell-Yan process, qq — > 7* — > from the NLO 

calculation for each mass bin in the rapidity region 
1.2 < \y\ < 2.2. C stands for the integrated luminosity, 
AMB/(oTot • £ mb C )i where Wmb stands for the number 
of sampled minimum-bias (MB) events and a^ ot repre- 
sents the total inelastic BBC MB cross section, 42.2 mb 
(2260 mb) for p+p (d+Au) collisions. eg/| c and eg^ c 
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FIG. 2: (color online) The invariant mass distributions between 5 GeV/c 2 and 16 GeV/c 2 are shown for p+p collisions (a,b) 
and in d+Au collisions for the South arm (c,d) and the North arm (e,f). In the mass distributions, (a), (c), and (e) show the 
mass distribution of unlike-sign foreground pairs, the mixed event pairs as combinatorial background, and the subtraction of 
background pairs from foreground pairs, (b), (d), and (f) show the combinatorial background subtracted signal overlaid with 
the correlated backgrounds and T(lS'+25'+35). The shaded bands around the curves represent the uncertainties from the 
fitting or calculations of the renormalization and factorization scales, Q/2 < n < 2Q. 
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FIG. 3: (color online) Differential cross sections for the Drell- 
Yan process, qq — > 7* — >• /i + /i~ are drawnfa) for a PYTHIA 
calculation and for NLO calculations [It], Gil in the rapidity 
region 1.2 < \y\ < 2.2. The ratios of the NLO calculations 
over that from pythia 6.4 (if-factor=1.4) are also shown (b). 



are the BBC trigger efficiencies for MB events and Drell- 
Yan events, respectively. C — £MB C / e DY C 1S a correction 
factor for the relative BBC efficiencies of minimum bias 
compared to hard processes containing a Drell-Yan pair. 
Its value is determined using a Glauber model and a sim- 
ulation of the BBC, and is 0.69 (0.89) for p+p (ri+Au) 
collisions. Ae^y represents the product of the detec- 
tor acceptance and efficiency, including the effect of the 
Level-1 trigger. Finally, N^y/Am is the yield of dimuon 
pairs from the Drell-Yan process for each mass bin. 

The detailed procedure to estimate the Drell-Yan 
yields for p+p and ri+Au, using Eq. ([2]) is as fol- 
lows. First, we generate the correct number of Drell-Yan 
events, which we estimate by multiplying the differential 
cross section by the accumulated luminosity for each in- 
variant mass bin considering BBC efficiencies. This cor- 
responds to dauy/dm ■ C ■ £q Y C m ^q. ©• For example, 
C is 22.5 pb _1 for the forward-rapidity p+p data. Af- 
ter event generation, to account for the acceptance times 
efficiency, Ae^y, the generated Drell-Yan events from 
the luminosity-weighted NLO calculation are then run 
through the PHENIX GEANT3 simulation and are recon- 
structed in the same way as real data. In the simulation, 



hit positions in all the muon detectors are registered and 
are reconstructed, including the effects of disabled HV 
channels and detector efficiencies. The resulting simu- 
lated counts in mass bins, N^y / Am, are then fit with an 
exponential function. This function describes the simu- 
lated distribution very well with a fit quality of \ 2 per 
degree-of-freedom (x 2 /dof) of 34.9/36 to 38.9/36. The 
shape and yield of this function are then fixed and used 
in the fits to the data, and represent the contribution of 
Drell-Yan in the fit function, Eq. ([5]). 

Although the NLO cross sections for ri+Au collisions 
already include nuclear corrections, they are still per 
nucleon-nucleon collision, and need to be scaled up by 
the number of binary collisions, iV co n. Eq. (JT]) shows the 
relation between the cross section for p+p and that for 
ri+Au collisions, which is derived from Eq. ([5]). 



RdAu — 1 



r1N dA 
DY 



n /dm 



(N co n}dNS Y /dm , 



(5) 



T riAu 

7 Tot 
PP 
7 Tot 



(dN^/dm) 



dAu 
7 Tot 



j riAu 
°DY 



I dm 



(JVcou) (dN^/dm) 



-pp 

7 Tot 



(Aeon) da^ Y jdm 



da^jdm = (JVcou) da™ Y /dm 



T dAu 
T Tot 
PP 
7 Tot 



(6) 
(7) 



where dN, 



riAu(pp) 
DY 



/ dm is the invariant yield for the Drell- 
Yan process in ri+Au (p+p) collisions and ct t ^ ( u ^ pp '' is the 
total inelastic cross section for ri+Au (p+p) collisions. 



In the expansion of Eq. (JB]), dN- t 



<iAu(pp; 
DY 



I dm 



riAu(pp) 

7 Tot 



is the differential cross section for the Drell-Yan process, 
I dm. (iV C oii) is the mean number of binary colli- 



der 



riAu(pp) 

DY 



sions and is calculated using a Glauber model and a simu- 
lation of the BBC. (iV C oii) is 7.6 ± 0.4 for inclusive ri + Au 
collisions. Eq. ([7]) is used for the Drell-Yan estimates in 
ri+Au collisions and da^y/^m is considered as nuclear- 
effect-corrected cross sections per nucleon-nucleon before 
scaling up. 



2. Correlations of open heavy-flavor pairs 

Several measurements of open bottom and charm 
cross sections have been made by PHENIX. A re- 
cent single-electron measurement of heavy-quark pro- 
duction at midrapidity obtained rj C e=551 ± 57 (stat) 
± 195 (syst) fib 23] for the total charm cross sec- 
tion. A dielectron measurement of the continuum 
charm pairs showed a total cross section of er C c=544 
± 39 (stat) ± 142 (syst) ± 200 (model) /Lib [Hj]. 
A perturbative-quantum-chromodynamics fixed-order- 
next-to-leading-log calculation [24| predicts a cross sec- 
tion of 256 ]ti4g /^b, which is within experimental and 
theoretical uncertainties with these measurements. 



8 



TABLE I: Simulation parameter settings for open 
beauty(charm) production. We used pythia 6.4 with 
the CTEQ5L parton distribution functions [2lj]. The pythia 
tunes are from a PHENIX dilepton mass spectra study [2^] 



Name of parameter Setting 

Bottom (Charm) Quark production on 

Bottom (Charm) Quark mass 4.1 (1.25) GeV/c 2 

k T 1.5 GeV/c 

if-factor 3.4 

Q 2 4 GeV 2 



Finally, a random angular correlation of the two open 
bottoms that form a pair is considered as an extreme 
case for the bottom correlation since NLO effects or in- 
teractions with other particles could alter the muon di- 
rections and destroy the angular correlation of the two 
heavy quarks. The pr spectra of the single muons from 
open heavy-flavor decay are sampled and the azimuthal 
correlation angle <f> of the decay muons is randomized 
for each muon and then pairs are formed, effectively de- 
stroying the angular correlation. The resulting difference 
between pythia estimation and random correlation is 
assigned as an additional systematic uncertainty. 



Existing bottom cross section measurements from 
PHENIX also agree within their uncertainties. An 
electron-hadron charge correlation measurement showed 
a total bottom cross section of a b ^=3.2 t\ \ (stat) 
-13 (syst) fjb [251 ] and a continuum mass distribution 
study obtained 0^=3.9 ± 2.5 (stat) ±1 (syst) /xb 22]. 
Meanwhile the calculation [24[ predicts ovr=1.87 
-o 67 ^ Hi consistent with the measurements. 

To obtain an estimate of the mass shape and to gener- 
ate simulated charm- and bottom-pair background events 
(see Tabic Q] for the simulation settings), we use the 
PYTHIA 6.4 tune, the same as that used for the PHENIX 
dilepton mass spectrum study [22| . The generated events 
are run through the PHENIX GEANT3 simulation to ac- 
count for the detector acceptance at forward-rapidity and 
are reconstructed by using identical code to that used to 
reconstruct the data, with the detector efficiencies in- 
cluded. Before reconstructing the simulated events, they 
are embedded into real events in order to match and eval- 
uate the effects of the multiplicity that exists for the data. 

The resulting mass spectrum is fit by an exponential 
function with (x 2 / dof ) of 10.0/12 to 11.8/12, which is 
then used to represent the open heavy-flavor component 
in the fits to the data, as in Eq. ©. The shape is fixed 
and the normalization, or yield, is allowed to vary in the 
fits to the data. The data fit-values obtained are within 
a half sigma in the experimental uncertainties of the pre- 
viously measured charm and bottom cross sections [221 ] 
(see Section Hi CI for further detail on the fits). For the 
shape of the fit function, we assign a systematic uncer- 
tainty by varying the slopes by ±10% from the nominal 
values obtained from the simulation. 

The relative ratio of bottom and charm production is 
fixed according to the measured production cross sec- 
tions obtained in the PHENIX dilepton mass spectrum 
study [22| with cr C c=544 ± 39 (stat) ± 142 (syst) ± 200 
(model) fib and a b ^=3.9 ± 2.5 (stat) (syst) fib. Since 
these have large measurement uncertainties, we assign a 
systematic uncertainty for the relative ratio of bottom 
and charm production cross sections by varying this ra- 
tio by ±100% from the nominal value; however, this does 
not result in a significant difference for the T yield be- 
cause the charm contribution is negligible in the T-mass 
region. 



3. T estimation 

Since the PHENIX muon-arm mass resolution is not 
good enough to resolve the three states of the T as shown 
in Fig. [4j we use results from two experiments at differ- 
ent collision energies and at different colliding systems 
to obtain an estimate of the relative ratio of the three 
T states, for the purpose of getting a distribution of line 
shape versus mass. The first is E605 [26], a p + A fixed 
target experiment at y/s NN =38.8 GeV and the second is 
CDF [27], a collider experiment with p + p collisions at 
y/s =1.8 TeV. These two experiments measured almost 
the same balance of the three T states although their 
energies and collision types are quite different (Table HI)). 



TABLE II: Relative strength of the three T states at CDF j27| 
and from FNAL E605 0. 



T states 


Mass 
(GeV/c 2 ) 


Branching ratio 
T ->■ n+fjT 


P + P [27] 
1.8 TeV 


p + A [26J 
38.8 GeV 


T(1S) 


9.46 


2.48% 


73% 


72% 


T(2S) 


10.02 


1.93% 


17% 


19% 


T(3S) 


10.36 


2.29% 


10% 


9% 



As the energy of the measurement reported here is be- 
tween those of the other measurements, we assume here 
that the composition of the three T states follows the 
ratio from CDF and assign a systematic uncertainty by 
varying the relative strength of the T(15 l ) over 0.73±0.10, 
with the fractions for the T(2S I ) and T(35) changing ac- 
cordingly. This uncertainty also accounts for the possi- 
bility that the 2S and 3S states might be suppressed more 
strongly than the IS in ci±Au minimum-bias collisions, 
since it allows for a 30% reduction in the 2S and a 50% 
reduction in the 3S. 

T simulations are performed for the three T states in 
order to estimate the effective T mass resolution and 
peak position for the real detector as well as to deter- 
mine the acceptance-times-efficiency correction. To ob- 
tain these estimates, we generate the three T states with 
pythia 6.4 and then process the generated T events 
through the PHENIX GEANT3 simulation to make events 
with hits in the detectors. These simulated events are 
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FIG. 4: (color online) T(1S+2S'+3S) are generated using 
pythia 6.4 and run through the PHENIX geant3 detector 
simulation. Each T has mass resolution of about 0.6 GeV/c 2 . 
The line shape used for the fit function, Eq. (J5J is composed 
of three Gaussians for the three T states. 
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tical code to that used to reconstruct the data - including 
resolution smearing effects. The sum of three Gaussian 
functions is fit to the T-mass distribution. The results 
from the fit provide an estimate of the widths and the 
means of masses for the three T states. The resulting 
shape, as shown in Fig. |U is then implemented for the 
fit function, Eq. © used to extract the yields from the 
data. 

Unlike the J/ip, where the mass resolution is predom- 
inantly determined by effects from multiple scattering 
in the absorber preceding the muon-tracking volume, the 
higher momentum muons from the T experience less mul- 
tiple scattering and less bending in the magnetic field; so 
the position resolution in the tracking volume becomes 
more important. To evaluate this, an additional system- 
atic is obtained by allowing the mass resolutions of the 
three states to vary by ±100 MeV/c 2 from their nominal, 
simulation determined, values. 

With the simulated T(lS l +25'-|-3S') events, we also cal- 
culate the acceptance times efficiency (Aer) by dividing 
the reconstructed T yields by the pythia generated T 
yields. Figure [3] shows Aex as a function of rapidity. Ae~f 
of the summed T(lS , +25'-|-3S') and of each T states sep- 
arately are quite similar to each other as shown. In this 
analysis, two inclusive rapidity bins are used, one for the 
positive rapidity and one for the negative rapidity. The 
values in d+Au collisions from 2008 are 0.0950 ± 0.0004 
and 0.0980 ± 0.0004 for positive and negative rapidity, 
respectively. For p+p collisions, in the rapidity same or- 
der, from 2006 they are 0.1132 ± 0.0007 and 0.1096 ± 
0.0007; and from 2009, they are 0.1164 ± 0.0007 and 
0.0907 ± 0.0007. 



FIG. 5: (color online) Acceptance x Efficiency of each T state 
and for the sum of the three states, T(lS+25'+3S'). The 
values for each state and for the sum of the three states are 
very similar. The T(IS+2S+3S) is comprised of the T(16 r ), 
T(25), and T(3S), following the abundances of 73% : 17% : 
10% from the CDF experiment See also Table |TI] 



then embedded into real events to reflect the multiplicity 
environment of the data and are reconstructed with idcn- 



C. Data evaluation and T extraction 

We extract Ts from the data using the estimated cor- 
related backgrounds and the Ts as described in Sec- 
tions III B 11 III B 21 and III B 31 The fit function used for 
the mass distribution is shown in Eq. ([5]). In addition to 
the T(15 l +25'-|-35') signal, the function includes contri- 
butions from the Drell-Yan process and from correlated 
open bottom/charm pairs. 



J 



F(m) = p exp(pim + p2m 2 + pzm 3 ') + pi[{l — p e )exp(p 5 m) + p 6 exp(p7m)] + 
Ps[(l -Pn -pi 4 )exp(-0.5A 2 ) + Pll cxp(-0.5B 2 ) +pi 4 cxp(-0.5C 2 )] 

I 



(8) 



where m is the invariant mass of the dimuon, A — (m — 
Pd)/pi0, B=(m -Pi2)/pi3, and C=(m - p 15 )/p le . The 
parameters po to p% are for the NLO Drell-Yan process 
and are fixed by the NLO calculations and the PHENIX 
GEANT3 simulation, as discussed in Section Hi B II 



The parameters p4 to pr are for the contribution of 
the open bottom and open charm correlated pairs. The 
relative ratio of bottom to charm yields, which is repre- 
sented by pe, is fixed from the PHENIX dilcpton mass 
spectra study [22j, and the shape is determined from 
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TABLE III: Parameter settings for the fits to the data and extraction of the Ts [Eq. JBJ], See text for the details. 



Parameter Fitting Parameter Setting 

po Yield of Drell-Yan process Fixed by NLO calculation 

Pi Slope of Drell-Yan process Fixed by NLO calculation 

P2 Slope of Drell-Yan process Fixed by NLO calculation 

P3 Slope of Drell-Yan process Fixed by NLO calculation 

P4 Yield of Charm/Beauty correlations Set free 

P5 Slope of Beauty correlation Fixed by pythia/geant simulation 
ps Relative ratio of Charm/Beauty correlations Fixed PHENIX dilepton measurement [22] 

P7 Slope of Charm correlation Fixed by pythia/geant simulation 

p 8 Yield of T(1S+2S+3S) Set free 

P9 Mean value of T(IS') Fixed by pythia/geant simulation 

pio Resolution of T(1S) Fixed by pythia/geant simulation 

pn Relative ratio of T(2S) Fixed by CDF experiment [27] 

P12 Mean value of T(2S) Fixed by pythia/geant simulation 

P13 Resolution of T(2S) Fixed by pythia/geant simulation 

P14 Relative ratio of Y(3S) Fixed by CDF experiment [27] 

pis Mean value of T(3S) Fixed by pythia/geant simulation 

pie Resolution of T(3S) Fixed by pythia/geant simulation 



TABLE IV: Systematic uncertainties for each source and for each collision type (see the text for details). Type A represents 
a point-to-point uncorrelated systematic uncertainty. Type B represents a common systematic uncertainty between points at 
different rapidity. Type C is a global uncertainty. 



Systematic uncertainty sources p+p South P+P North d+Au South d+Au North Type 

(backward rapidity) (forward rapidity) (backward rapidity) (forward rapidity) 



Relative ratio of T(1S+2S+3S) 


2.6 % 


0.8 % 
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0.3 % 


0.1 % 
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T mass resolution 


6.4 % 


6.7 % 
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8.3 
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NLO DY model 


7.7 % 


6.9 % 


4.0 
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NLO DY renormalization/factorization 


0.9 % 


0.9 % 


7.7 


% 
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Open Bottom random correlation 


7.2 % 


4.0 % 
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% 


6.3 


% 
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Combinatorial background normalization 


0.6 % 


0.7 % 
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% 
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Combinatorial background estimation 


1.7 % 
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0.0 
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MuID efficiency 
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BBC efficiency 
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the PYTHIA and the PHENIX GEANT3 simulation, as de- 
scribed in Section III B 21 The total yield from correlated 
bottom and charm, p 4 , is allowed to vary in the fits to 
the data. The resulting contributions from correlated 
bottom and charm are then checked against those from 
the PHENIX dilepton measurements, which have bot- 
tom and charm cross sections of cr C c=544 ± 39 (stat) 
± 142 (syst) ± 200 (model) [ib and 0^=3.9 ± 2.5 (stat) 
(syst) /xb [22j]. For this check, we integrated our fit- 
ted correlated bottom and charm over the mass range 5 



GeV/c 2 to 16 GeV/c 2 , and then added the contribution 
from unmeasured regions assuming the mass shape from 
the NLO calculation. This estimate of the contribution 
of bottom and charm is within a half sigma in the exper- 
imental uncertainties of the nominal cross sections. For 
d+Au collisions, this estimate is still within a half sigma 
when we scale the nominal cross section by the number 
of participants (2x 197) assuming no nuclear modification 
effects on the production. 
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of T, is allowed to vary in the fits to the data in order 
to extract the T signal. 

The data are fit using a log-likelihood fitting method 
that adds the normalized combinatorial background to 
both the mass distribution and the fit function. This 
has the advantage that empty bins in the mass distri- 
bution, which result from statistical fluctuations of the 
background above the signal size and otherwise produce 
negative counts, are accounted for properly. The fitting 
quality is very good with x 2 /dof of 9.0/16 and 6.4/16 
for the negative and positive rapidities of p+p collisions, 
respectively, and 14.6/16 and 9.5/16 for the negative and 
positive rapidities of d+Au collisions, respectively. Sys- 
tematic uncertainties from the NLO calculation, the as- 
sumed cross sections, and detector performances are ex- 
plained in Sections III B 1 1 Hi B 21 and III B 3l and are sum- 
marized in Table IIVI Figure [6] (a), (b), and (c) show 
the T(1S , +2S'+3>S') mass distribution after subtracting 
off all correlated backgrounds, for p+p and d+Au col- 
lision systems. We checked above the T-mass region, 
> 11.5 GeV/c 2 , for fitting reliability; the integral of the 
high-mass region is within the systematical uncertainties, 
which are drawn as shaded bands. 
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FIG. 6: (color online) The mass distributions are drawn 
after subtraction of all the correlated backgrounds and 
Drell-Yan process. The fitted peak curves represent only 
T(1S'+2S+3S'). The shaded bands (green vertical shading) 
around the fitted peak curves represent fitting uncertainties 
and those around zero yield (blue slanted shading) represent 
systematic uncertainties from the different assumptions for 
NLO Drell-Yan models, pythia and random correlations of 
charm/bottom. 



The parameters p% to p\§ are for the contribution of 
the T(15+2S+35). A, B, and C represent the means 
of masses and widths of the three T states, as estimated 
using the PHENIX GEANT3 simulation package and fixed 
for data - as described in Section III B 31 The total yield 



TABLE V: T invariant yields and cross sections of p+p and 
d+Au data sets are shown. The first uncertainty shown is the 
statistical, and the second is the systematical. 



Data set N^l 



(xlO 1 



£W-diVr/dy 
(xl(T 10 ) 



B„ 



do"T /dy 



p+p South 5.1 3.8 ± 1.1 ± 0.5 16.0 ± 4.6 ± 2.2 pb 

p+p North 5.2 3.6 ± 1.0 ± 0.4 15.4 ± 4.3 ± 1.8 pb 

d+Au South 1.3 17.8 ± 5.5 ± 2.9 4.0 ± 1.2 ± 0.6 nb 

d+Au North 1.3 25.2 ± 5.6 ± 3.3 5.7 ± 1.3 ± 0.7 nb 



The invariant yields of T(15+2S'+3>S') for each rapid- 
ity bin are calculated as, 



N r ■ C 



dN r _ 
dy Nmb ■ Aer ■ Ay ' 



(9) 



where the notation is the same as for Eq. (|4]), except T 
is used instead of the Drell-Yan process. Table [V] shows 
calculated yields and cross sections for both arms and for 
both collision types. 

Figure [7] shows the invariant T yields for p+p and 
d+Au collisions. The nuclear modification factor, RdAu, 
can be obtained from the invariant yields. Eq. ([5]) shows 
the relation between RdAu, the invariant yield, and N co u. 
The scale factor, N co \\, makes RdAu one if the T yield for 
d+Au collisions is equal to the T yield for p+p collisions 
times the number of binary collisions in ri+Au collisions, 
i.e. RdAu =1 if there are no nuclear modification effects. 
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Rapidity 

FIG. 7: (color online) Invariant yields, B^^dNr /dy, for p+p 
and d+Au collisions are shown as a function of rapidity. The 
solid error bars represent the statistical uncertainties. The 
boxes represent the systematic uncertainties. The global sys- 
tematic uncertainties are quoted as text at the bottom. 
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FIG. 8: (color online) Nuclear modification factors, RdAu, are 
shown as a function of rapidity. For comparison, the upper 
panel (a) shows RdAu for the J/tp [|| as well as for the T. The 
solid error bars represent the statistical uncertainties and the 
boxes represent the systematic uncertainties. The global sys- 
tematic uncertainty is quoted as text at the bottom. The 
lower panel (b) shows theoretical predictions of nuclear mod- 
ification based on a NLO EPS09 combined with a breakup 
cross section, with <T{, r =0 to 8 mb in 2 mb steps from top to 
bottom. See text for the details. 



As seen in Fig. [8j at forward rapidity T production 
shows no significant suppression with an RdAu of 0.91 
± 0.33 (stat) ± 0.16 (syst); while at backward rapidity 
the suppression of the T is approximately one sigma (of 
the experimental uncertainty) below one with an RdAu 
of 0.62 ± 0.26 (stat) ± 0.13 (syst). Figure|U(a) shows a 
comparison to previous results from PHENIX for RdAu 
of the J/ip- The J ftp results show a larger suppression at 
forward than at backward rapidity, a trend that cannot 
be confirmed or denied for the T given the large uncer- 
tainties of the measurements presented here. 

A NLO calculation with EPS09 shadowing and a 
breakup cross section [lOj predicts modest suppression 
at backward rapidity, but no shadowing at forward ra- 
pidity; although there could be suppression by a breakup 
cross section, as seen in the red lines in Fig. [5] (b). The 
rapidity dependence of this NLO calculation appears to 
be consistent with the trend between our backward- and 
forward-rapidity measurements. At both backward and 
forward rapidities, the large uncertainties of the measure- 
ments do not give a significant constraint on the breakup 
cross section within the context of the NLO models. It 
will also be of interest to compare other models that in- 
clude the effects of initial-state parton energy loss or of 
gluon saturation to this T data. 

We can parameterize the nuclear dependence of T pro- 
duction as cr^ +Au = o-P +p x (2A Au ) a for d+Au colli- 
sions, where Aa u represents the number of nucleons in 
the gold nucleus. As for RdAu, if there are no nuclear ef- 
fects then a would be one. Previously, E772, which was 
at y/ s NN =38.8 GeV, showed a large decrease in a at xf 

< 0. The PHENIX backward-rapidity covers -0.42 < x F 

< -0.14, where x F — x\ — X2 and x\ is the momentum 
fraction of the gluon in dcuteron. The backward-rapidity 
((x F ) ~ -0.2, (x 2 ) ~ 2 x 10" 1 ) PHENIX measure- 
ments obtain ar £(lS+2S+ 3S) = 0-^25 ± 0.070 (stat) 
± 0.035 (syst) and at forward rapidity ((xp) ~ 0.2, (a^) 
- 1 x 10~ 2 ) aT(lS , +2S , +3S') = °- 990 ± °- 060 ( stat ) 
± 0.029 (syst). 

Figure [5] shows a versus xf and versus x 2 from the 
E772 data and from our data. The suppression levels 
of T(lS'+25+35') in PHENIX are consistent with those 
from E772 within uncertainties. 

For our c?+Au measurements, we can also calculate 
the ratio of the T yield, dN^^ u /dy, between back- 
ward and forward rapidities as a test of the nPDF. This 

ratio, d^ +Au /dy|-2.2< y <-i.2 / dN^ +Au ldy\ x . 2<y<2 . 2l 
shows some suppression at backward rapidity relative to 
forward rapidity, with a value of 0.71 ± 0.27 (stat); but 
the effect is not very significant due to the large uncer- 
tainty. 



IV. SUMMARY & CONCLUSION 

In summary, we have presented the first yields, cross 
sections, and nuclear dependences for T production in 



13 



0.95 
0.9 
0.85 



0.65 
0.6 



1.05 
1 

0.95 
0.9 

0.85 
0.8 

0.75 
0.7 

0.65 



(a) 



■ T(1S),^S^ = 38.8GeV(E772) 

YPS+SS),^!^ 38.8 GeV(E772) Z 
• T(1St2S+3S),^S^ = 200GeV(PHENIXj; 



X 



(b) 



r 1 



■ T(1 S), = 38.8 GeV(E772) 

T(2S+3S), ^S,^ = 38.8 GeV(E772) 
• T(1S+2S+3S),\(S BN = 200GeV(PHENIX) 

10' 



FIG. 9: (color online) a versus xf and xi. The square and 
diamond shaped points are from the E772 experiment where 
an 800 GeV proton beam collides with fixed targets of 2 H, 
C, Ca, Fe and W, corresponding to \Js NN =38.8 GeV. The 
round points are from this analysis. The solid error bars rep- 
resent the statistical uncertainties and the boxes represent the 
systematic uncertainties. 



y/s NN =200 GeV d+Au and p+p collisions for two rapid- 
ity bins. At backward rapidity, T(1S+2S+3S) yields are 
measured to be suppressed by approximately one sigma 
of the experimental uncertainty below one. The rapid- 
ity dependence of the observed T suppression at forward 
and backward rapidities are compatible with lower en- 
ergy results and a NLO theoretical calculation. Com- 
parison to the theoretical calculation for a model that 
includes EPS09 shadowing and a breakup cross section 
does not result in any definitive constraint on the breakup 



cross section given the large experimental uncertainties. 
Future comparisons to gluon saturation models and to 
models including initial-state energy loss would also be 
of interest. 
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